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ABSTRACT 


A form of Maxwell' s moment method is used to obtain constitutive rela- 

2 

tions for stress and heat flux for a gas passing through a shock wave. The c 

X 

moment of the Boltzmann equation together with several physical approximations 
are used to obtain a nonlinear expression relating stress and deformation rate, 
with Mach number as a parameter. It is shown that the constitutive relations 
for both stress and heat flux reduce to the Navier-Stokes relations for the case 
M = 1 and |r/p| « 1, and differ from these relations by as much as a factor 
of two for |t/p| =1. Theoretical density profiles for shock-wave Mach numbers 
between 2. 8 and 8 are compared with experimental data from shock waves in 
argon and the agreement is shown to be excellent, thus supporting the use of the 
proposed constitutive relations. The physical significance of the approximations 
used and areas for further improvement in the theory are discussed. 
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I. INTRODUCTION 


An important step in the application of the moment method employed by 

Mott-Smith to analyze shock structure, and which was later extended by Lees 

to treat problems in rarefied gas dynamics, is the selection of an approximate 

1 2 

functional form for the distribution function ’ . The functional form selected 

must satisfy physical requirements imposed by boundaries and possess sufficient 

flexibility, in the form of free parameters, to satisfy a desired number of moments 

of the Boltzmann equation. The power of the method resides in the fact that for 

situations involving extreme translational nonequilibrium one can often construct 

a distribution function, relying mainly on intuition and a physical understanding 

of the problem, which is a better approximation to the actual flow than one 

obtained through a more rational approach such as the Chapman- Enskog procedure. 

Mott-Smith' s bimodal distribution and Lees' two-sided Maxwellian are well- 

known examples of distribution functions constructed on this basis. 

% 

A useful interpretation for explaining the effectiveness of this particular 
approach to Maxwell’s moment method is that it provides, in effect, a mechanism 
for closing a set of moments of the Boltzmann equation by expressing, the unknown 
higher order moments in terms of the "known" lower order moments through the 

3 

assumed form of the distribution function. For example, the moment <e > 

* A 

4 

(or <c x >) is the unknown moment in the Mott-Smith solution to the shock- 

structure problem and the bimodal distribution allows one, in principle, to express 
3 

<c > in terms of the lower order moments, i. e. , gas dynamic variables, 
x 


1 



Since the specific information needed to close a set of moments of the 

Boltzmann equation is one or more particular higher order moments of the 

distribution function, this raises the following question: Is it easier to guess the 

3 

moment needed, e. g. < c x >» or to guess the distribution f and then compute the 
moment needed? For problems such as Couette flow or heat conduction in a 
rarefied gas, the latter is probably most likely true since solid surfaces play 
such important roles, and it seems natural that processes occurring at a surface 
should enter directly into the assumed form for f. For the shock-structure 
problem, there is reason to believe the former is true because of the absence of 

3 

solid boundaries and because of the fact that one of the unknown moments <c > 

x 

is directly related to physical quantities in the flow; and, as will be seen below, 
one can therefore estimate this moment rather well without having to select an 
approximate f. 

The aim of this paper is therefore to present an analysis of the shock- 
structure problem with a basically different emphasis on the Maxwell moment 
method than that employed by Mott— Smith. In fact, the aim is to show that for 
steady, one-dimensional flow one can close the set of moments of the Boltzmann 
equation, at the "Mott-Smith" level of approximation, by approximating the 
unknown moments directly rather than approximating f. The primary reason 
for emphasizing this approach is that one arrives at new constitutive laws for 
stress and heat flux, which replace the Navier*- Stokes relations, and consequently 
all of the advantages that accrue from the continuum description of a gas are 
retained. In addition, it becomes possible to see directly the physical significance 


2 



of the solution and specifically where the analysis could be improved, in con- 
trast to previous solutions where the outcome of the analysis is simply a shock- 
wave profile for specific assumptions. 
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II. MOMENT EQUATIONS 


For a steady, one-dimensional situation, the "first" five moments of 
the Boltzmann equation are given by 


d? pu) - 0 


(la) 


d 2 

-£(p+pu -T) = 0 


(lb) 


~j“[P u (h +|u 2 ) +(q -tu)J = 0 


(lc) 


d?< C x> = mA[C x J 


(Id) 


d 2 2 2 

d?< C x C > = mAtC x C J ’ 


(le) 


where p,p,h,u,r and q are the standard gas dynamic variables, m is the 

2 2 

mass of the molecule, A[c x ] represents the c^ moment of the collision 
integral, and c x is the molecular speed in the x-direction (the direction of 
flow). The terms in Eqs. (Id) and (le) arising from the collision integral can 

3 

be expressed in terms of t and q for Maxwell molecules , i. e. , 


mA[c 2 ] = 

1 x J |1 


(2a) 


2 2p 2 

mA[c x C ] = -^TU-^q), 


(2b) 


where the constants describing the intermolecular potential have been replaced 
by the coefficient of viscosity p, as defined by the Chapman-Enskog procedure. 
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This simplification is frequently used in the approximate theories of rarefied 
gas dynamics. 

In order to show where the present work could be generalized in the 
discussion below, by removing the restriction to Maxwell molecules, it is 
convenient to define here the quantity r* by the relation 

mA[cJl = Jr*. (3) 

Use of Eq. (3) will therefore preserve the algebraic form obtained for Maxwell 
molecules since t* = t in that case. In general, one would have to obtain the 
relation between t* and r (and possibly other variables) for a particular inter- 
molecular potential. 

Equations (1) and (2) together with the relations for a thermally and 
calorically perfect gas 


p = pRT 

(4a) 

h = c T 
P 

(4b) 


do not form a complete set of equations since there are more unknowns (9) than 

3 2 2 

equations (7); that is, the moments <c > and <c c > are the extra unknowns. 

If the approach used by Mott-Smith were pursued at this point, then the extra 
unknowns would be evaluated from an assumed distribution function. Rather than 
follow this approach, let us first consider the information contained in Eqs. (Id) 
and (le). 

Using the relations for Maxwell molecules (2a) and (2b), we have from 
(Id) and (le) 
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(5a) 


p d 3 

iZ T ' S’* V 


2D 2 . d 2 2 

7*™ - V° = 5? <c x c 


Since t is determined by (5a) therefore q is essentially determined by (5b). If 
3 2 2 

the moments <c > and <c c > could somehow be expressed as functions of the 
xx 

gas dynamic variables, it is clear that Eqs. (5) would yield constitutive laws for 
stress and heat flux. It is therefore desirable to preserve the identities of (5a) 
and (5b) during all algebraic manipulations in order to identify the resulting 
constitutive relations, and in addition, to remind one that all approximations 
introduced will affect only these relations. Once these relations are established, 
it is a simple procedure to compute shock-wave profiles with the complete set of 
equations and thereby check the approximations used by comparing the results 
with experimental data. 

Since Eq. (5b) is an equation for the heat flux q, it would seem that a 

far better approximation for q than that obtained by guessing f and computing 
2 2 

<c c > would be to assume the total enthalpy is conserved within a shock-wave 

X 

profile and use the energy equation (lc) to obtain 

q=Tu (6) 

4 

Equation (6) has been shown by Liepmann, Narasimha, and Chahine to be an 

5 

excellent approximation for the Navier-Stokes equations, and Sherman and Talbot 
have shown the approximation is very good for several continuum theories as well 
as for the Mott-Smith solution. Here we will assume Eq. (6) holds in general and 


rely on the comparisons given below to justify this selection. 
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Using the energy equation consistent with (6) and the notation defined by 
(3), the set of equations given by (1) becomes 


d 

d? pU) = ° 

(7a) 

d 2 

^(p+pu -T) = 0 

(7b) 

d . 12 

d? h + i U > = 0 

(7c) 

T* J“p<C 3 > 

p dx x 

(7d) 

q = tu. 

(7e) 


Displayed in the above form, we see that the constitutive laws for stress and heat 
flux must arise from Eqs. (7d) and (7e) respectively. In fact, one should expect 
to recover the Navier-Stokes relations 


o 4 du 

= ys 


(8a) 


o 

q 



(8b) 


from these two equations as the strength of the shock wave approaches zero, i. e. , 

as the Mach number approaches unity. The nature of this limit process is not 

obvious when inspecting Eq. (7d); however, the expression can be cast in a form 

for which the limit process becomes evident immediately. This can be seen if 
3 3 

we express <c > in terms of a central moment <V >, where V is the thermal 

X XX 

speed given by 
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V = c - u. 

X X 


( 9 ) 


If we then let M ->-1 and f where f^ is the Maxwellian distribution, 

s 

3 (0) 3 

<V > becomes zero since r 1 is an even function of V and <V > is an odd 
x xx 

moment. Therefore the remaining terms in Eq. (7d) must reduce exactly to (8a). 

3 

A. Expression for < c x > 

3 3 

The relation between <c > and <V > can be found directly from the 

X X 

definition (9), i. e. , 


3 

c 

x 



+ U) 


3 


and 


3 , r 3 . 2 3 

<c > = <V > + 3u<V > + u , 
xx x 

since by definition <V > = 0. Recalling that the x-component of the pressure 

X 

tensor is defined to be 


p xx = P< V = P - T - 


the above relation becomes 


3 3 3 

p<c > = p<V >+3u(p-f)+pu . (10) 

X X 

3 

Equation (10) shows that <c > can be split into a sum of two parts, a central 

X 

3 

moment <V x > and terms involving gas dynamic quantities. 

Differentiating Eq. (10) with respect to x, expanding the right hand 
side, and using the continuity and momentum equations, (7a) and (7b), one obtains 
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d 3 d 3 2,du 

S p< V = 3? < V + [3<p - T) - pu te 

2 5 n 

Substituting this expression into (7d) and using the relation a = - - for the 

3 p 

isentropic speed of sound in a monatomic gas, we have 



where M is the local Mach number. The only difference between (7d) and (11) 

3 3 

is that <c > has been replaced by the central moment <V >, otherwise both 
x x 

equations are exact. 

The limit value for x* for a weak shock wave can now be found from 
Eq. (11). When the shock Mach number is very close to one, the local Mach 
number M is very nearly unity throughout the shock wave, jr/pj is much less 

3 

than one, and <V > is nearly zero. Therefore the right hand side of Eq. (11) 

reduces to which is precisely the Navier-Stokes expression. For a 

general molecular force law, it is not clear that t* would also reduce to t for 

M =1, but at least for Maxwell molecules it is certain that Eq. (11) yields the 
s 

Navier-Stokes relation (8a). 

This result immediately suggests that Eq. (11) could be quite useful for 
constructing an approximate expression for x* for larger Mach numbers. The 

3 

assumption <V > = 0 gives a very good result for M near one; thus any 

X S 

3 

reasonable approximation to <V > should improve the expression, at least for 

X 

a certain Mach number range. Here we can also see an advantage in using 
Eq. (11) over Eq. (7d) if we were to apply Mott-Smith' s approach at this point 
by guessing a suitable f. That is, the entire expression for x* in Eq. (7d) 
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would be approximate whereas only a part of the expression for t* in Eq, (11) 

would be approximate, and it would only be a minor contribution if M is near 

s 

one. In a certain sense one can view the correction as arising from a perturba- 
tion to the shock Mach number from unity. 

A further inspection of Eq. (11) shows that it, unfortunately, is not well 
behaved for all M g , since the term in square brackets becomes zero in the up- 
stream part of the shock wave for 5 = #75 = 1.34. Therefore for strong 

3 

shock waves the contribution to t* due to <V x > becomes very important to the 
total expression and it no longer can be considered as arising from a small 
perturbation. 

It is of interest that the critical Mach number predicted by Eq. (11), 

3 

when the term involving <V > is ignored, is close to the critical Mach 

g 

numbers for the various continuum theories; e. g. , Sherman and Talbot report 

that solutions of Grad' s thirteen-moment equations are limited to the range 

1 SM < 1. 65 and solutions to Burnett' s equations are limited to the range 
s 

ISM < 2. 1. In view of the fact that difficulty has been encountered with many 
s 

schemes, for even moderate Mach numbers, it can be seen that an approximation 

3 

for <V > must be chosen with care, 
x 

3 

B. Approximation for <V x > 

3 

Before proceeding with approximate schemes for evaluating <V >, the 

X 

3 

physical interpretation of p<V > should be stressed* Using the definition of the 

x 

x-component of the heat -flux vector, we have 
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( 12 ) 


, = |p<v x v 2 >. 

Obviously for a one-dimensional flow, such as a shock wave, only the x-component 
of the heat-flux vector is nonzero. However, this single component q can be 
split into two parts, if we define q" and q A by the relations 

q " = |p <V x> ( 13a ) 

qX = 2^ <V X ^ V y + V z )># ( 13b ) 

We can then speak of the parallel and perpendicular parts of the heat flux q; 
thus our problem requires an approximation for the parallel part of the heat 
flux as shown by (11) and (13a). 

Since q, q", and q x have identical dimensions, it is convenient to define 
the nondimensional quantity 

Q = q "/q, (14) 

\ 

which exhibits the ratio of the parallel part to the total heat flux q. If Q is 

3 

assumed given, then <V > is obtained from Eqs. (13a) and (14), i. e. , 

p<V^> = 2qQ. (15) 

The reason for introducing the ratio Q becomes apparent once one 
evaluates the quantity for the Chapman-Enskog distribution f^. For this case 
we have 

Q = 3/5, (16) 

i. e. , Q is strictly a constant everywhere in a flow field. The Maxwellian 
distribution f^ can be viewed as giving the "zeroth" approximation to <V^> 
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which is zero; and here r ^ yields what can be considered to be the "first" 
approximation, i. e. , 

P< v x> = (17) 

dT 

Since q in a certain sense varies like — , it is easy to see that the 

2 

second term in (11) is essentially — while the first term is thus it makes 

dx 

sense that the second term should vanish first for weak shock waves as discussed 
above. 

It is a straightforward procedure to substitute (17) into (11), use approxima- 
tion (6), and obtain a new expression relating stress and deformation rate for a 
gas passing through a shock wave. The algebraic effort is reduced considerably 
if the following substitution is carried out first 

1 dq Id t du [ u dr 

pdx pdx pdx p dx ’ 

d T 

where — is obtained from the momentum equation (7b). This leads to the 
useful result 


m 2 -w p & 


(18) 


However, it should be remembered that (18) is only an approximate relation 
since (6) was used in its derivation. 

g 

Using Eq. (18) together with the Chapman- Enskog value for p<V > given 
by (17) and substituting into (11), the "first" approximation for t* is found to be 
given by 
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( 19 ) 


t* = [f(l - 1 ) - ”M 2 jj±™. 
5' p 15 dx 


Here again 9 we see that for M g ~ 1 and |r/pj « 1 the limiting value of t* 
is However, for this approximation the critical Mach number is given by 

O UX 


M = ^27/7 1 = 1.96, 


which is an improvement over the previous value, but it is still not very exciting. 

The simplest generalization to the Chapman-Enskog value is to assume 
Q is still a constant but that it takes on a different value a, the value being 
determined by the vanishing of the critical Mach number. Following the same 
substitution as above, we now have 


T * = [(! -t/p) +|^~M 2 ](3 - 2o!)p~. 

2 

For a = 5/6 the coefficient of the M term becomes zero and the expression 
no longer exhibits a critical Mach number. Our new approximation is then 

Q = 5/6, (20) 

and we obtain a new law relating stress and deformation rate given by 

* ,, / ,4 du 

T»=<l-T/p)j|X S . 

A more symmetric form of the equation is given by 

T ° =t*/(1 -t/p), (21) 

where t° is the Navier-Stokes value (8a). 

For Maxwell molecules r* = r and Eq. (21) together with (7c) and (7e) 
yield the extremely simple constitutive relations 
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T = r/(l - t/p) 


q = q/(l - T/p), 


(22a) 


(22b) 


where in the second equation we have assumed (J-c^/k = 3/4, consistent with 
approximation (6), in order that the Navier-Stokes relations are fully recovered 


for small values of 


M- 


A plot of Eq. (22a) is presented in Fig. 1, where the dimensionless strain 
rate T°/p is plotted against the dimensionless stress r/p. Since t/p is always 
a negative number in a shock wave, only one quadrant of the coordinate system 
is shown in the figure. The upper limit on the value of |x/pj is well defined as 

4 

discussed by Liepmann, Narasimha and Chahine ; i. e. , as a result of the 
definition of x/p kinetic theory limits its value to -2, and manipulation of the 
momentum equation shows that its value is limited to -3/2 which occurs for an 
infinitely strong shock wave. The figure therefore shows that the values of 
-T°/p predicted by Eq. (22a) are always less than the Navier-Stokes values and 
the ratio to the Navier-Stokes value reaches a minimum of 2/5. Consequently, 
shock wave thickness as predicted by (22a) should be greater than that predicted 
by the Navier-Stokes relations. 

Although the constitutive relations (22a) and (22b) are remarkably simple, 
which suggests a certain validity in itself, it is appropriate to question the 
approximation (20) in view of the disagreement with the Chapman-Enskog value, 

Q = 3/5. However, since the Chapman-Enskog value is in effect the result of a 
special limit process, it' s quite conceivable that other values could be obtained 
for other conditions. 
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For very large shock Mach numbers, the flow at the upstream end of a 

shock wave where significant molecular collisions first begin to occur is akin 

to a molecular beam, because of the developing difference between the parallel 

and perpendicular components of the thermal motion of the molecules. If we 

2 

consider an actual molecular beam as an example, for which <V > / 0 and 

x 

2 

<Vy> = 0, then definitions (13a), (13b) and (14) show that Q takes on its 
maximum value for this case and its value is given by 

Q = q"/(q" +q i ) = 1, (24) 

where q A = 0. Likewise the definition shows that the absolute minimum value of 
Q would be zero. Therefore one would expect Q to assume values between these 
two extremes in an actual flow field. 

In fact, it would even be difficult to imagine that Q could be less than 1/3 
in a shock wave, since this would mean that the transport of energy associated 
with each transverse degree of freedom exceeds the transport of energy associated 
with the parallel degree of freedom. If the value 1/3 is truly approached some- 
where in a shock wave, it seems reasonable to expect to find it at the subsonic 
end of a strong shock wave where the Mach number is low and the distribution 
function loses its preferential direction. 4 

Likewise, the upper limit Q = 1 is probably not attained in a shock wave 
since the upstream portion is never precisely a molecular beam regardless of 

how large M becomes; and It is quite possible that the requirement that the 
s 

critical Mach number vanish in Eq, (11), which led to the result Q = 5/6, in 

fact selects the maximum permissible value of Q in a shock wave. 

If we now assume that Q lies in the range 

15 



1/3 < Q < 5/6, 


(25) 


with large values on the supersonic side and small values on the subsonic side of 
a shock wave, it becomes clear that the Chapman-Enskog value Q = 3/5 could 
very likely be the proper value for the portion of the flow where the local Mach 
number is unity. 

The next level of approximation to (15) is therefore to establish a functional 
form for Q beyond merely a constant. To do this one must identify the important 
independent variables entering the expression for Q. We have already suggested 
that the Mach number is important is arriving at (25). The constitutive relations 
show that t/p is a very important quantity in a shock wave, so it is likely that it 
too is an important variable in expressing Q, If we chose to add to this list the 
second Chapman-Enskog parameter q/pa, we find however that it is not a new 
variable in a shock wave, since it is a combination of M and t/p, i. e. , 

q/pa = Tu/pa = (t/p)M. (26) 

Since Eq. (26) lists the most important variables in a shock wave, it is reasonable 
to assume that a suitable approximation would consist of the following expression 
Q = Q(M,t/ P ), (27) 

with (25) setting the upper and lower limits on Q for all M and t/p. 

Several important details concerning the functional form of Q could now 
be established by substituting (15) into (11), using the assumed form (27), and 
studying the various limits of the resulting expression. With this information, 
a particular algebraic expression could be chosen for (27) (this would be equivalent 
to choosing f); and thus an improved approximation for t* would then be obtained. 
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The important feature of this approximation would be that T°/p would depend on 
both t/p and M (or q/pa) rather than merely on t/p as given by (22a) ; this 
would be the first level of approximation at which the constitutive law for stress 
would also depend on the heat flux q» as one would expect in a nonequilibrium 
situation. Since the algebraic details for this case can become involved, it is 
appropriate to defer this discussion and first obtain a prediction of shock-wave 
profile based on the simpler constitutive relations (22a) and (22b). 
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III. SHOCK- WAVE STRUCTURE, Q = 5/6 


The constitutive relations (22a) and (22b) were obtained using approxima- 
tion (6), assuming Maxwell molecules, a Prandtl number of 3/4, and the approxima- 
tion Q = 5/6. The analysis to this point can be viewed as an effort to derive ap- 
proximate constitutive relations for stress and heat flux for a gas passing through 
a shock wave, not to solve the full set of Eqs. (1). However, we can now return 
to the set given by (1) and replace Eqs. (Id) and (le) with the approximate rela- 
tions (22a) and (22b). In fact, since it is reasonable to assume that a material 
description of a gas cannot be very sensitive to the type of gas being considered 
(for example, the Navier-Stokes relations apply equally well to both monatomic 
and diatomic gases, as well to liquids), one can postulate that the above constitu- 
tive relations hold in general and solve the set for the appropriate viscosity law 
p(T), ratio of specific heats y, and Prandtl number Pr needed to describe the gas. 

Using this point of view, the appropriate set of equations to solve for shock 
structure is given by 


pu = m 

(28a) 

m(u - u x ) +p - p x =r 

(28b) 

m(h +|u 2 - H 1 ) = (tu - q) 

(28c) 

T° =t/(1 -t/p) 

(28d) 

q° = q/(l - t/p) 

(28e) 

o 4 du 

T = ~ 11— 

3dx 

(28f) 
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(28g) 



to which one would add (4a) and (4b) for a thermally and calorically perfect gas. 

In the above set, the first three equations are integral forms of Eqs. (la) to (lc), 
subscript 1 is used to indicate an upstream value, and p and k are assumed to 
be known functions of temperature. 

As a result of the fact that the above expressions for t and q yield the 
following relation 
o,o , 

t /q = T/q, (29) 

we will find that one can use exactly the same method to solve the above set as that 

g 

discussed by Gillarg and Paolucci for the solution of the Navier-Stokes equations. 
In fact, Eq. (29) guarantees the same behavior at the upstream and downstream 
singular points for the solution in the T, u plane as that found by Gillarg and 
Paolucci for the Navier-Stokes equations. 

Before proceeding with the analysis, it is quite useful to introduce an 
algebraic tracer e as a coefficient before the t/p term in both Eqs. (28d) and 
(28e) and then use the following scheme: If we set e = 0, the above set of 
equations become identical to the Navier-Stokes equations; and if we set e = 1, 
we recover the set of Eqs. (28). 

Using the algebraic tracer e , we have from Eq. (28d) 

T°/p = (t/p)/(1 - £t/p) = (Tu/pu)/(l - CTu/pu). 

Multiplying the momentum Eq. (28b) by u, to obtain an expression for tu, and 
defining the dimensionless variables T = T/T^ and u = u/u^, the above equation 
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can be written in the following form 


t°/ Pi = T0/u(T - e0), (30) 

where 

e = [T - 5(1 +yM 2 (i - u))]. 

S 

Simplification of the equation for the two cases e = 0 and e = 1 is obviously 
straightforward. Likewise, from Eq. (28d), we have 

q°/pu = (q/pu)/(l - e-ru/pu). 

Using the energy Eq. (28c) to obtain q and the momentum Eq. (28b) to obtain 
tu, the above equation can be written in the dimensionless form 

q /p lUl = t[ 0 +^~(i + V M s d ' u > " T)]/(T " e0) * (31) 

If we divide Eq. (31) by Eq. (30) to obtain the differential equation for the solution 
curve in the T, u plane, we see that the algebraic tracer e cancels out, i. e. , 
we obtain 

= -|Pr(y - 1)M 2 [0 +7771 +I f i,M s ( 1 - S 2 )) " T]u/0. (32) 

Therefore the integral curve in the T, u plane is exactly the same as that for the 
Navier-Stokes solution. However, the velocity gradient and the temperature 
gradient are quite different for the two solutions as shown by Eqs. (30) and (31) 
for the two cases e = 0 and e = 1. 

Once Eq. (32) has been integrated to obtain the universal curve T = T(u) 

for a particular combination of M , y and Pr, Eq. (30) can be evaluated to 

s 
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obtain the dimensionless velocity gradient t 0 / p M , for either the present 

i s 

solution (e = 1) or the Navier-Stokes solution (e = 0). Figure 2 presents a 
comparison of the two solutions for several different shock Mach numbers and for 
the case y = 5/3, Pr = 2/3. The figure shows that the two solutions are nearly 
identical for M = 1. 5. However, they are quite different for the stronger shock 
waves, with the present solution having a significantly smaller value for the 
maximum velocity gradient. It can be seen from the figure that most of the 
change occurs on the supersonic side of the shock wave (u * 1) with very little 
change occurring on the subsonic side. 

In order to generate a shock wave profile from the data in Fig. 2, it is 
necessary to specify the temperature dependence of the coefficient of viscosity P- 

U> 

Here we assume a power law dependence for the viscosity, p ~ T , with 
co = 0. 75 (to represent argon as discussed below). Figures 3 and 4 give a 
comparison of several density profiles for the two solutions. The spatial 
variable x has been nondimensionalized by dividing by the upstream mean-free 

i6 rr ^1 

path length \ = -r-*/- . The two profiles were arbitrarily matched at the 

Jl O If 

subsonic side because of the suggestive appearance of the curves in Fig. 2. 

The figures show that the two solutions yield nearly identical results for 
M g = 1.5 and the difference in the two profiles becomes greater as the shock 
Mach number increases. The increase in shock-wave thickness over the Navier- 
Stokes value shown by the present solution is definitely in the proper direction 
for agreement with experimental results. However, we also see from Fig. 4 
that considerable asymmetry develops in the predicted density profile for the 
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larger values of M g ; and this is not in agreement with experimental results. 
If a detailed comparison is made with experimental density profiles, one finds 
that the present solution lies too close to the Navier-Stokes solution at the 
subsonic end of the shock wave, i. e. , we must conclude that the constitutive 
relations require further corrections for |t/p| small. 
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IV. EXTENSION, Q = Q(M) 


The first step beyond the approximation Q = 5/6 must include the most 
important functional dependence in the representation for Q. The discussion 
leading to Eq. (27) suggests that M and t/p are important variables in repre- 
senting Q, but no judgment was made as to which variable is the most important. 
By assuming initially that both must be included, we can arrive at an equation, 

Eq. (36), which will show that the functional dependence on M is far more 
essential than the dependence on t/p. 

Assuming the representation Q = Q(M, t/p), the x-derivative of the 
product qQ is given by 

(33) 

where Q M and are partial derivatives with respect to M ''and t/p re- 

spectively. Using the definitions for M and the isentropic speed of sound a, 
we have 

dM _ d u _ 1 du u_ d_ j yp 

dx ~ dx'a ~ a dx 2 dx V p ’ 

a 

which by virtue of the continuity Eq. (28a) allows one to write 


dM 1 du Id. 
dx - a dx 2p dx PU ” 

Multiplying the momentum Eq. (28b) by u and using approximation (6) together 
with Eq. (18), the above expression becomes 
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Similarly for the x-derivative of t/p, we can write 


ir /p = l? q/pu) = - < T/p >|? pu >i- 

Following the same procedure as above for determining the x-derivative of pu, 
we have finally 

d, l r/ll -2 2.2 / ,du 

jpr/p.jKM -l)+(jM +DT/P15J. (35) 

Using the derivatives given by (34) and (35) together with the expression 
derived in (18) for Eq. (33) can be written as 

p i? )Q = I q<m2 ■ 1 +t/p> + vy +i > Mx/p ^/ p [(m2 ■ i) 

+ (|m 2 + x )t /piJ|, 

where the approximation q = ru has again been used. With the aid of definition 
(15) and the above equation, the expression for t* given by Eq. (11) yields the 
result 

4t*/t° = [(9 - 6Q) - (5 - 6Q)M 2 ] - [(9 - 6Q) - 2Q M (M 2 +3)M 

(36) 

- 2Q r/p (M 2 - l)]r/p +[2Q r/p (2M 2 +3)](r/p) 2 . 

Since |t/p| s 3/2 for all shock Mach numbers and since in general 
t* ». t, therefore t*/t° must be bounded; and this fact allows several important 


conclusions to be drawn concerning the functional form of Q. If it were assumed 



that Q is a function of t/p alone, then the last two bracketed terms in Eq. (36) 
would increase with M 2 and t*/t° would not be bounded. Therefore Q cannot be 
a function of t/p alone. However, if it were assumed that Q is a function of 

O 

M alone, then it would only be necessary to require to vary like M for 
large M and t*/t° would be bounded. Therefore the dependence on M is 
necessary in the functional representation for Q, as long as Q is not a constant. 

Proceeding on the basis that Q = Q(M) is the appropriate generalization 
beyond the previous approximation Q = 5/6, Eq. (36) simplifies to read 

4t*/t° = [(9 - 6Q) - (5 - 6Q)M 2 ] - [(9 - 6Q) - 2Q M (M 2 + 3)M]t/ P . (37) 

Returning to the requirement that t*/t° be bounded, one concludes from (37) that 
Q(M) must exhibit the following behavior for large M 

M 2 t5 - 891 = »! <38a) 


lim 
M -oo 



(38b) 


where c and c are constants. If we further assume c and c are nonzero, 

12 1 a 

the functional form of Q(M) for large M can be established in greater detail. 

In particular, Eq. (38a) shows that the difference [5/6 - Q] must vary like M . 
If in addition we agree that Q(l) must be the Chapman-Enskog value 3/5 and 
that Q(M) is a monotonic function, we can write 

Q(M) = 15M 2 /(7 +18M 2 ) (39) 

as a reasonable approximation for Q(M) satisfying the cited requirements. 

It is obvious that Eq. (39) satisfies Eqs. (38a) and (38b) since 
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[5/6 - Q] ss 35/6(7 +18M 2 ), 

-3 

and therefore Q^. ~ M which is required by (38b), In arriving at Eq. (39) it 
was assumed that the simplest algebraic expression satisfying the above condi- 
tions would be adequate. 

Figure 5 presents a plot of the function Q(M) defined by Eq. (39). The 
curve was not drawn for M < 1/2 since the minimum Mach number behind a 
shock wave in a monatomic gas is approximately 1/2. Although the proposed 
function has the correct behavior for large M (in the sense that c and c are 

I Cl 

nonzero), it is not obvious that it has the correct behavior near M = 1, which 
unfortunately cannot be determined with the information at hand. Of particular 
concern is the large variation exhibited by the proposed function about the 
Chapman-Enskog value, which possibly is a result of the simple algebraic 
expression chosen. 

The validity of the approximation could be checked quite conclusively by 
including in Fig. 5 data obtained from numerical solutions of the Boltzmann 
equation or any suitable approximation thereto. The numerical data, however, 
would have to be obtained for a wide range of shock Mach numbers and from solu- 
tions for which the shock wave is accurately described in the wings of the profile 
(where the Mach number reaches either bounding value) in order to adequately 
check the proposed function. Numerical data in a preliminary form was very 
kindly made available to the authors recently by S. M. Yen who together with B. 

L. Hicks has used the Monte Carlo method to evaluate the Boltzmann collision 
integral for hard sphere collisions. The general trend of this preliminary data 
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has provided considerable encouragement for the use of Eq. (39), and in par- 
ticular it has helped to dispel certain misgivings which were originally held 
about the large variation in the function near M = 1 from the Chapman-Enskog 
value. (Because of the preliminary nature of the data, it was not felt that it would 
be appropriate to include these data in Fig. 5. ) 

Turning our attention back to Eq. (37) and using Eq. (39), we again find a 
relatively simple constitutive relation for stress, i. e. , 

T ° = t*/(F - Gt/p), (40) 

where F and G are functions of M given by 

F = (63 +37M 2 )/4(7 +18M 2 ) 

G = 3(147 +126M 2 +292M 4 )/4(7 +18M 2 ) 2 . 

Constitutive relations for t and q can again be obtained by assuming the gas is 
composed of Maxwell molecules and by following the same algebraic manipulation 
leading to Eqs. (22a) and (22b). It is of interest, and of importance to the 
development below, that the only change that occurs in Eqs. (22) is the quantity 
(1-t/p) is replaced by (F - Gt/p). Since F and G are bounded by the values 

1.57 > F >0.357 

1.42 >G >0.676 

for the Mach number range 0. 5 £ M £ 00 , one can see the modification is not 
one of a fundamental nature. 

The equation relating stress and deformation rate for Maxwell molecules, 

obtained from Eq. (40), is plotted in Fig. 6 for fixed values of M. Here we see 

the Navier-Stokes relation is recovered only for the case M = 1 and |r/p j 
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small. The fact that the Navier-Stokes relation is not recovered far upstream 
and far downstream from the shock wave where |t/p| is small and M / 1 
comes as a surprise from both a physical and a mathematical point of view. This 
difference could, in part, be attributed to the approximation used for Q and the 
fact that the r/p dependence was ignored. However, the difference appears in 
precisely the correct direction to remove the asymmetry from the shock-wave profile 
that was noticed in Figs. 3 and 4 with the simpler theory. This fact in itself is a 
strong argument in favor of employing the constitutive relation for stress as 
displayed in Fig. 6 and checking the accuracy with which the new theory can 
predict shock-wave profile. 

Since Eq. (29) still holds and Eq. (40) introduces M as the only new 
variable, the same method can be used to compute shock-wave structure with the 
new approximation as was discussed in the above section for the simpler theory; 

Eq. (32) is of course exactly the same and Eqs. (30) and (31) must be modified 
in a trivial manner to include the functions F and G appearing in Eq. (40). 

Figures 7 and 8 present the results of several such calculations for the 
density profile through a shock wave for the case y = 5/3, Pr = 2/3, and 
co = 0. 75. A suitable value of co for low temperature argon is 0. 816 which has 
been used by a number of investigators since the work of Gillarg and Paolucci. 

If shock-wave profiles are to be computed for M g s 10, then a value is need 
for a much greater temperature range. For the temperatures T 2 10, 000°K, 
a suitable "average" value appears to be co = 0. 75, according to the calculations 

7 

for argon by Amdur and Mason and the experimental data obtained by Camac 
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8 9 10 

and Feinberg , Collins and Mennard , and Willeke . The dashed 

curves appearing in the two figures are plots of experimental data for argon 
obtained by B. Schmidt 11 using an electron beam densitometer and the 17-inch- 
diameter shock tube at Caltech, The theoretical and experimental density pro- 
files were arbitrarily matched at the point (p - p )/(p - p ) = 0. 5. 

X A 1 

Two observations become immediately evident on comparing the theoretical 
and experimental curves in the figures. First, the theory predicts the experimental 
data rather well over the entire Mach number range from 2. 8 to 8; and second, the 
theoretical profiles no longer exhibit the asymmetry that was present in Figs. 3 
and 4. A close inspection of the figures shows that the theoretical profiles are 
slightly thinner than the experimental profiles in each case and the fit is slightly 
better for the larger Mach numbers. The value of oj used in the calculations 
could certainly be adjusted to improve the fit still further, especially for the 
lower Mach numbers, but in view of the uncertainty associated with the proper 
value of w for a wide temperature range in argon and in view of the further im- 
provements that could be incorporated in the present theory, it seems advisable to 
judge the theory on the basis of the present comparison alone. On this basis it is 
quite clear that the constitutive relation for stress exhibited in Fig. 6 is more 
nearly correct than that displayed in Fig. 1, and therefore any further improve- 
ment in the theory would probably yield results in closer agreement with Fig. 6. 

This again raises the very interesting question of why the Navier-Stokes relations 
are not recovered far ahead of the shock wave and far behind the shock wave where 

the flow is nearly in equilibrium, since most physical arguments would tend to 
favor the opposite conclusion. 
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V. DISCUSSION 


In view of the favorable comparisons exhibited by Figs. 7 and 8, further 
work to confirm or improve the present theory would be clearly justified on the 
grounds that it would contribute to the macroscopic interpretation of the physical 
processes taking place within a shock wave. For example, it would be extremely 
interesting, from the point of view of continuum mechanics, to check the constitu- 
tive relation for stress plotted in Fig. 6. This could be done with data obtained 

from numerical solutions of the Boltzmann equation or from calculations involving 

12 

collisions between a finite number of particles, such as the work by Bird. The 
calculations would have to be carried out for several different values of M 

s 

since data corresponding to a single shock- wave profile would appear as a single 

loop in Fig. 6, starting from the origin and ending at the origin, with the curve 

for M = 1 dividing the loop into two nearly-equal halves. Displayed in this form, 

the data should show a nonlinear dependence on t/p and a Mach-number dependence 

about the locus of points for M = 1. If Fig. 6 is approximately correct, the 

trajectory of each loop near the origin, for large values of M , would have a 

s 

larger slope for the supersonic branch than for the subsonic branch. In terms of 
the present interpretation of p(T), this implies that the Navier-Stokes relation 
for stress does not apply at the two ends of a shock wave. 

A possible resolution of the problem could come from a more careful 
treatment of Eq. (40), e.g. , for |r/p| « 1 we have 
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( 41 ) 


o 4 du 
= Ft = -F^T)-. 

The procedure has been to assume Maxwell molecules (i. e. , r* = t) and then to 
assume p.(T) does not have to be restricted to the temperature dependence T 
characterizing Maxwell molecules; obviously the two sides of Eq. (41) are not 
treated equally in such a procedure. An important extension to the present work 
would therefore be to relate t* to t (and possibly q) for a general inverse- 
power intermolecular potential so that compatible expressions for r* and p(T) 
could be used; and thus physical interpretation of (41) would become more meaning- 
ful. 

Figure 5 would likewise provide a physically interesting area for comparison 
with numerical solutions. Namely, the quantity Q expresses the ratio between 
1) the rate of flux of thermal energy contained in the degree of freedom which is 
parallel to the direction of propagation of a shock wave and 2) the total rate of 
flux of thermal energy. This information is of course useful in constructing a 
picture of how the distribution function changes from an upstream to a downstream 
position. In the case of Fig. 5. , it would be desirable to determine whether the 
representation Q = Q(M,t/p) would be necessary or whether the present approxima- 
tion Q = Q(M) is sufficient. 

A final check on the present theory would be to review the approximation 
q = tu, or equivalently H = constant, used in deriving the constitutive relations 
for stress and heat flux. The present approach is at least self consistent in that 
the maximum variation in the quantity (H - H^)/H^ for all of the profiles shown 
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in Figs. 7 and 8 is about 2%. However, this is not surprising since the integral 
curve T = T(u) is the same as that for the Navier- Stokes solution, and therefore 
the relation H = H(u) should also be the same. 

The constitutive relation for stress has been discussed in terms of the 
independent variables r/p and M rather than the equivalent combination of 
variables r/p and q/pa (see Eq. (26)), because of the fact that one can recover 
the case M = 1 so much more easily. In terms of a physical description of a 
gas, it would seem that r/p and q/pa would be a more meaningful combination, 
since the velocity of the fluid does not appear in either variable. This raises 
the interesting possibility that the constitutive relations written in this form 
could also be applied to nonsteady, one-dimensional flows, e. g. , to study the 
reflection of a shock wave from a wall, to analyze the nonsteady development of a 
shock front, or to study the attenuation of ultrasonic waves in a gas. 
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Fig. 1. Relation between nondimensional stress t/p and nondimensional 
strain-rate T°/p (for Q = 5/6) compared with the Navier-Stokes 
relation. 




3 


2 

-TtyPiMs 


1 


0 


Fig. 2. 



Nondimensional velocity gradient for the Navier-Stokes solution and 
for the present solution (for Q = 5/6). 
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Fig. 6. Constitutive relation for stress based on the approximation Q = Q(M), 
defined by Eq. (39). 
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